######################################################################
# Time-stamp: <05/08/11 16:11:59 ostolop>
######################################################################
package EP::EPC::Core::Clustering::exp;

use strict;

use EP::Config;

use EP::Common;
use EP::Common::Data;
use EP::Common::General;

use EP::Common::DBI::Dataset;

use File::Temp;
use File::Basename;

use EP::Common::Error ':try';
use EP::Rclient;

######################################################################
# process expression data clustering
######################################################################
sub process_clustering {
  my $self = shift;
  my ( $dataset, $query ) = @_;

  my $matrix = $dataset->data->raw_data;

  $_ = $query->param ( "clustering_action" );

  my $action = { action_id       => 0,
     visualizable    => 1,
     action_date     => gmtime,
     dataset_id      => $dataset->id,
         };

 SWITCH: {
    /hierarchical/ && do {
  $dataset = hierarchical ( $dataset, $query );
  $action->{component_id} = 3;
  $action->{action_metadata} = { distance           => $dataset->{distance},
               transpose_distance => $dataset->{transpose_distance},
               rows_clust   => $dataset->{clustering},
               cols_clust   => $dataset->{transpose_clustering},
               heatmap_data => $dataset->data->{matrix},
               rows => $dataset->data->{row},
               cols => $dataset->data->{column},
               distance_measure => $query->param ( "distance" ),
               algorithm => $query->param ( "algorithm" ),
               collapse_genes => $query->param ( "COLLAPSE_GENES" ) || "",
               collapse_genes_percent => $query->param ( "COLLAPSECRITERIA_GENES" ) || "",
               collapse_conditions => $query->param ( "COLLAPSE_CONDITIONS" ) || "",
               collapse_conditions_percent => $query->param ( "COLLAPSECRITERIA_CONDITIONS" ) || "",
#              layout => $query->param ( "layout" ),
               impute_row_avg => $query->param ( "impute_row_avg" ),
               fast_nr_dists => $query->param ( "fast_nr_dists" )
             };

  last SWITCH;
      };
    /kmeans/
      && do {
  $dataset = kmeans ( $dataset, $query );
  $action->{component_id} = 4;
  $action->{action_metadata} = { distance           => $dataset->{distance},
               clustering         => $dataset->{clustering},
               heatmap_data => $dataset->data->{matrix},
               rows => $dataset->data->{row},
               cols => $dataset->data->{column},
               K => $query->param ("K"),
               distance_measure => $query->param ( "distance" ),
               start => $query->param ( "start" ),
               impute_row_avg => $query->param ( "impute_row_avg" )
             };
  last SWITCH;
      };
    /kmedoids/
      && do {
  $dataset = kmedoids ( $dataset, $query );
  $action->{component_id} = 4;
  $action->{action_metadata} = { distance           => $dataset->{distance},
               clustering         => $dataset->{clustering},
               heatmap_data => $dataset->data->{matrix},
               rows => $dataset->data->{row},
               cols => $dataset->data->{column},
               K => $query->param ("K"),
               distance_measure => $query->param ( "distance" ),
               start => $query->param ( "start" ),
               impute_row_avg => $query->param ( "impute_row_avg" )
             };
  last SWITCH;
      };
    /similarity_search/
      && do {
  $dataset = similarity_search ( $dataset, $query );
  $action->{component_id} = 13;
  $action->{action_metadata} = { distance           => $dataset->{distance},
               clustering         => $dataset->{clustering},
               heatmap_data => $dataset->data->{matrix},
               rows => $dataset->data->{row},
               cols => $dataset->data->{column},
               K => $dataset->{K},
               N => $query->param ( "N" ),
               genes => $query->param ( "genes" ),
               distance_measure => $query->param ( "distance" ),
               impute_row_avg => $query->param ( "impute_row_avg" ),
               merge_sim_results => $query->param ( "merge_sim_results" )
             };
  last SWITCH;
      };
  }

  $dataset->{action} = $action;

  return 1;
}

sub hierarchical {
  my ($dataset, $query) = @_;
  my $cluster_what = $query->param('cluster_what');
  my $fast_nr_dists = $query->param('fast_nr_dists');

  $fast_nr_dists = 100 if not $fast_nr_dists;

  my $custom_distances_rows = $query->upload ( "custom_distances_rows" ) if $query->param ( "custom_distances_rows" );
  my $custom_distances_cols = $query->upload ( "custom_distances_cols" ) if $query->param ( "custom_distances_cols" );

  my $custom_distances_rows_fh;
  my $custom_distances_cols_fh;

  if ( $custom_distances_rows ) {
    $custom_distances_rows_fh = new File::Temp ( DIR => "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location, UNLINK => 0, SUFFIX => '.dist' );
    print $custom_distances_rows_fh $_ while <$custom_distances_rows>;
  } else {
    $custom_distances_rows_fh = 0;
  }

  if ( $custom_distances_cols ) {
    $custom_distances_cols_fh = new File::Temp ( DIR => "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location, UNLINK => 0, SUFFIX => '.dist' );
    print $custom_distances_cols_fh $_ while <$custom_distances_cols>;
  } else {
    $custom_distances_cols_fh = 0;
  }

  $fast_nr_dists = 100 if $custom_distances_rows or $custom_distances_cols;
  $fast_nr_dists = 100 if $query->param("distance") !~ /(eucl|manh)$/;

  if ( $query->param ( "distance" ) ne "none" ) {
      if ( $cluster_what != 1 ) {
    # rows or rows and columns
    $dataset = create_distance_matrix ( $dataset, $query->param ( "distance" ), $query->param ( "impute_row_avg" ), $custom_distances_rows_fh ) if $fast_nr_dists == 100;
    $dataset = create_hierarchical ( $dataset, $query->param ( "distance" ), $query->param ( "algorithm" ), $fast_nr_dists );
      }

      if ( $cluster_what != 0 ) {
    # columns or rows and columns
    $dataset = create_transpose_distance_matrix ( $dataset, $query->param ( "distance" ), $query->param ( "impute_row_avg" ), $custom_distances_cols_fh ) if $fast_nr_dists == 100;
    $dataset = create_transpose_hierarchical ( $dataset, $query->param ( "distance" ), $query->param ( "algorithm" ), $fast_nr_dists );
      }
  }

  return $dataset;
}

sub kmeans {
  my ($dataset, $query) = @_;

  $dataset = create_kmeans ( $dataset, $query->param ( "distance" ), $query->param ( "K" ), $query->param ( "start" ), $query->param ( "impute_row_avg" ) );

  return $dataset;
}

sub similarity_search {
  my ($dataset, $query) = @_;

  $dataset = create_similarity_search ( $dataset, $query->param ( "distance" ), get_gene_ids ( $dataset, $query->param( "genes" ), $query->param("rows_regexp") ), $query->param("N"), $query->param ( "impute_row_avg" ), $query->param ( "merge_sim_results" ) );

  return $dataset;
}

sub create_transpose_distance_matrix {
  return create_distance_matrix ( @_, "transpose" );
}

sub create_transpose_hierarchical {
  return create_hierarchical ( @_, "transpose" );
}

sub create_hierarchical {
  my ($dataset, $distance, $algorithm, $fast_nr_dists, $transpose) = @_;

  if ( $fast_nr_dists == 100 ) {
    my ($dist_file)  = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/$dataset->{distance}";
    my ($clust_file) = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".$dataset->{distance}.$algorithm.cluster";

    if ( defined $transpose ) {
      $dist_file  = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/$dataset->{transpose_distance}";
      $clust_file = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".$dataset->{transpose_distance}.$algorithm.transpose_cluster";
    }

    if ( -f $clust_file && ( -M $clust_file < -M $dist_file ) ) {
      # clustering is already performed
      warn "[CLUSTERING] Using previously calculated clustering results";
    } else {
      # in theory could contain multiple ones.... allows us to do parallel jobs
      my $jobs = { 0 => { program => $EP::Config::h->{HIERARCHICAL_CLUSTERING_CMD},
        param => "-m $algorithm -d $dist_file -o $clust_file",
        stdout => "$clust_file.comments",
        stderr => "$clust_file.comments"
      } };
      die "Failed to run the clustering algorithm!" if ( EP::Common::General::submit_jobs_and_wait ( $jobs ) ) != 1;
    }

    $dataset->{clustering}           = $dataset->filename . ".$dataset->{distance}.$algorithm.cluster"                     if not defined $transpose;
    $dataset->{transpose_clustering} = $dataset->filename . ".$dataset->{transpose_distance}.$algorithm.transpose_cluster" if     defined $transpose;
  } else {
    if ( $algorithm eq "ave" ) {
      $algorithm = "average";
    } elsif ( $algorithm eq "min" ) {
      $algorithm = "single";
    } elsif ( $algorithm eq "max" ) {
      $algorithm = "complete";
    } else {
      $algorithm = "average";
    }

    my $clust_file = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".$algorithm.cluster";
    $clust_file = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".transpose.$algorithm.cluster" if defined $transpose;

    unlink $clust_file if DEBUG;

    if ( -f $clust_file ) {
      # clustering is already performed
      warn "[CLUSTERING] Using previously calculated fast clustering results";
    } else {
      my $rows = $dataset->data->raw_data->dim(1);
      $rows = $dataset->data->raw_data->dim(0) if defined $transpose;

      warn "Rows: $rows, fast: $fast_nr_dists";
      my $nr_dists = int ( 0.01 * $fast_nr_dists * ( $rows * $rows / 2 ) ); # calculate about $fast_nr_dists percent of distances

      $nr_dists = 250000 if $nr_dists > 250000;

      # HACK -> quickly create a NetCDF to test things out
      my $dataset_nc = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".nc";
      $dataset_nc = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".transpose.nc" if defined $transpose;

      if ( not -f $dataset_nc ) {
        use PDL::NetCDF;
  my $ncobj = new PDL::NetCDF ( $dataset_nc );
  if ( defined $transpose ) {
    $ncobj->put ( 'matrix', ['arrays', 'genes'], $dataset->data->raw_data->setbadtoval(0)->transpose );
  } else {
    $ncobj->put ( 'matrix', ['genes','arrays'], $dataset->data->raw_data->setbadtoval(0) );
  }
  $ncobj->close;
      }

      # in theory could contain multiple ones.... allows us to do parallel jobs
      my $jobs = { 0 => { program => $EP::Config::h->{FAST_HIERARCHICAL_CLUSTERING_CMD},
        param => "data_in=$dataset_nc dist=$distance cluster=$algorithm pivstr=quantile small=$nr_dists random=$nr_dists clust_out=$clust_file",
        stdout => "$clust_file.comments",
        stderr => "$clust_file.comments"
      } };
      die "Failed to run the fast clustering algorithm!" if ( EP::Common::General::submit_jobs_and_wait ( $jobs ) ) != 1;
    }

    $dataset->{clustering}           = $dataset->filename . ".$algorithm.cluster"           if not defined $transpose;
    $dataset->{transpose_clustering} = $dataset->filename . ".transpose.$algorithm.cluster" if     defined $transpose;
  }

  return $dataset;
}

sub create_kmeans {
  my ($dataset, $distance, $K, $start, $impute_row_avg ) = @_;

  $dataset->{distance} = $distance;

  my ($bin_file)  = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".bin";
  my ($clust_file) = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".$dataset->{distance}.K=$K.kmeans.start=$start.cluster";

  $bin_file = nan_impute_row_avg ( $bin_file ) if $impute_row_avg;

  if ( -f $clust_file ) {
    # clustering is already performed
    warn "[CLUSTERING] Using previously calculated clustering results in $clust_file";
  } else {
    # in theory could contain multiple ones.... allows us to do parallel jobs
    my $jobs = { 0 => { program => $EP::Config::h->{KMEANS_CLUSTERING_CMD},
      param => "IN=$bin_file OUT=$clust_file K=$K dist=$distance -start start=$start",
      stdout => "$clust_file.comments",
      stderr => "$clust_file.comments"
          } };
    die "Failed to run the clustering algorithm!" if ( EP::Common::General::submit_jobs_and_wait ( $jobs ) ) != 1;
  }

  $dataset->{clustering}           = $dataset->filename . ".$dataset->{distance}.K=$K.kmeans.start=$start.cluster";
  return $dataset;
}

sub create_kmedoids {
  my ($dataset, $distance, $K, $start, $impute_row_avg ) = @_;

  $dataset->{distance} = $distance;

  my ($bin_file)  = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".bin";
  my ($clust_file) = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".$dataset->{distance}.K=$K.kmedoids.start=$start.cluster";

  $bin_file = nan_impute_row_avg ( $bin_file ) if $impute_row_avg;

  if ( -f $clust_file ) {
    # clustering is already performed
    warn "[CLUSTERING] Using previously calculated clustering results in $clust_file";
  } else {
    # in theory could contain multiple ones.... allows us to do parallel jobs
    my $jobs = { 0 => { program => $EP::Config::h->{KMEDOIDS_CLUSTERING_CMD},
      param => "IN=$bin_file OUT=$clust_file K=$K dist=$distance -start start=$start",
      stdout => "$clust_file.comments",
      stderr => "$clust_file.comments"
          } };
    die "Failed to run the clustering algorithm!" if ( EP::Common::General::submit_jobs_and_wait ( $jobs ) ) != 1;
  }

  $dataset->{clustering}           = $dataset->filename . ".$dataset->{distance}.K=$K.kmedoids.start=$start.cluster";
  return $dataset;
}

sub create_similarity_search {
  my ($dataset, $distance, $gene_list, $N, $impute_row_avg, $merge_sim_results ) = @_;

  $dataset->{distance} = $distance;

  # need to add 1 to each gene id for Jaak's K-means business (1-based)
  $_++ foreach @$gene_list;

  my $start = "list:" . join ( ",", @$gene_list );
  my $maxnr = $N + 1;   # algorithm returns input gene as one of the found genes
  my $K = scalar @$gene_list;

  my ($bin_file)  = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".bin";
  my ($clust_file) = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".$distance.K=$K.N=$N.$$.cluster";

  $bin_file = nan_impute_row_avg ( $bin_file ) if $impute_row_avg;

  # in theory could contain multiple ones.... allows us to do parallel jobs
  my $jobs = { 0 => { program => $EP::Config::h->{KMEANS_CLUSTERING_CMD},
          param => "IN=$bin_file OUT=$clust_file K=$K dist=$distance -start start=$start -select and maxnr=$maxnr",
          stdout => "$clust_file.comments",
          stderr => "$clust_file.comments"
        } };
  die "Failed to run the clustering algorithm!" if ( EP::Common::General::submit_jobs_and_wait ( $jobs ) ) != 1;

  $dataset->{clustering}           = $dataset->filename . ".$distance.K=$K.N=$N.$$.cluster";
  $dataset->{K} = $K;
  return $dataset;
}

sub get_gene_ids {
  my ($dataset, $inputGenes, $regexp) = @_;

  my $dataset_folder_path = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $dataset->folder->location;

  my @inGenes = split /\r\n|\r|\n/, $inputGenes;
  warn "Looking for @inGenes" if DEBUG;

  my @gene_ids;
  my $rows = $dataset->data->rows;

  foreach my $row ( keys %{$rows} ) {
    my $full_annot = $rows->{$row}->{full_annot};
    my @annot_names = $full_annot->Keys;

    foreach my $ingene ( @inGenes ) {
      if ( $regexp ) {
  if ( grep { $_ =~ /$ingene/} $rows->{$row}->{full_annot}->Values ) {
    push @gene_ids, $row;
    warn "Found gene $ingene at row $row!" if DEBUG;
  }
      } else {
  if ( grep { $_ eq $ingene } $rows->{$row}->{full_annot}->Values ) {
    push @gene_ids, $row;
    warn "Found gene $ingene at row $row!" if DEBUG;
  }
      }
    }
  }

  die "Couldn't find $inputGenes in the dataset!" if @gene_ids == 0;
  return \@gene_ids;
}

sub kmedoids {
  my ($dataset, $query) = @_;

  my $custom_distances_fh = $query->upload ( "custom_distances" ) if $query->param ( "custom_distances" );

  $dataset = create_distance_matrix ( $dataset, $query->param ( "distance" ), $query->param ( "impute_row_avg" ), $custom_distances_fh );
  $dataset = create_kmedoids ( $dataset, $query->param ( "distance" ), $query->param ( "K" ), $query->param ( "start" ) );

  return $dataset;
}

sub nan_impute_row_avg {
  my ($bin_file) = @_;
  my $nanimp_bin_file = $bin_file . ".nanimp";

  return $nanimp_bin_file if -s $nanimp_bin_file;

  EP::Rclient->R(<<__RCODE__);
app_root = "$EP::Config::h->{ROOT}";
source("$EP::Config::h->{ROOT}/cgi-lib/R/ep.io.R");
source("$EP::Config::h->{ROOT}/cgi-lib/R/ep.impute_data.R");

M = ep.readBin("$bin_file");
O = ep.nanimprowavg(M);
ep.writeBin( O, "$nanimp_bin_file" );
__RCODE__

  die "R did not generate imputed data!" if not -e $nanimp_bin_file;
  return $nanimp_bin_file;
}

sub create_distance_matrix {
  my ( $dataset, $distance, $impute_row_avg, $custom_distances_fh, $transpose ) = @_;

  if ( "$custom_distances_fh" ne 0 ) {
    warn "Using custom distance matrix $custom_distances_fh!";

    $dataset->{distance} = basename ( "$custom_distances_fh" )           if $transpose ne "transpose";
    $dataset->{transpose_distance} = basename ( "$custom_distances_fh" ) if $transpose eq "transpose";

    return $dataset if $custom_distances_fh;
  }

  my ($bin_file)  = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".bin";
  my ($dist_file) = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".$distance.dist";

  if ( $transpose eq "transpose" ) {
    $bin_file  = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".transpose.bin";
    $dist_file = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/" . $dataset->filename . ".$distance.transpose.dist";
  }

  $bin_file = nan_impute_row_avg ( $bin_file ) if $impute_row_avg;

  my ($rows)    = $dataset->{row_count};
  my ($columns) = $dataset->{column_count};

  if ( -f $dist_file && ( -M $dist_file < -M $bin_file ) ) {
    # distance matrix is ready for use
    warn "[CLUSTERING] Using previously calculated distances";
  } else {
    # in theory could contain multiple ones.... allows us to do parallel jobs
    my $jobs = { 0 => { program => $EP::Config::h->{DISTANCES_CMD},
      param   => "dist=$distance IN=$bin_file OUT=$dist_file",
      stdout  => "$dist_file.comments",
      stderr  => "$dist_file.comments"
          } };

    die "Failed to run the distance matrix calculation!" if ( EP::Common::General::submit_jobs_and_wait ( $jobs ) ) != 1;
  }

  $dataset->{distance}           = $dataset->filename . ".$distance.dist"           if $transpose ne "transpose";
  $dataset->{transpose_distance} = $dataset->filename . ".$distance.transpose.dist" if $transpose eq "transpose";

  warn "GOT $dataset->{transpose_distance}";

  return $dataset;
}

1;
